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Abstract A non-equilibrium steady state can be characterized by a nonzero but stationary flux driven by a static 


external force. 


Under a weak external force, the drift velocity is difficult to detect because the drift motion is feeble 


and submerged in the intense thermal diffusion. In this article, we employ an accurate method in molecular dynamics 
simulation to determine the drift velocity of a particle driven by a weak external force in a one-dimensional periodic 
potential. With the calculated drift velocity, we found that the mobility and diffusion of the particle obey the Einstein 
relation, whereas their temperature dependences deviate from the Arrhenius law. A microscopic hopping mechanism 
was proposed to explain the non-Arrhenius behavior. Moreover, the position distribution of the particle in the potential 
well was found to deviate from the Boltzmann equation in a non-equilibrium steady state. The non-Boltzmann behavior 
may be attributed to the thermostat which introduces an effective “viscous” drag opposite to the drift direction of the 


particle. 
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1 Introduction 

Non-equilibrium phenomena are ubiquitous in the 
universe. In particular, many machines, including 
biomolecules in living cells, often function in a non- 
equilibrium steady state (NESS).!'-4] When a system is 
driven by a static “force” (such as electric field, shear, 
chemical-potential gradient, etc.) and dissipates at the 
same time to reach a stationary state, it works in an 
NESS,?:4] which is quite common and pivotal in many ap- 
plications, such as propellants, batteries, sensors, electro- 
spray, and electrophoresis systems. Therefore, the prop- 
erties and processes of systems working in an NESS are 
of significant fundamental and practical interests. One 
unique feature of an NESS is the nonzero but station- 
ary flux driven by an external force, which characterizes 
the active transport process in the system. The response 
of a system to a strong external force is noticeable and 
the unidirectional drift velocity is easy to detect in both 
experiments®—®] and simulations.2~-!8] However, in most 
applications, the system usually works in a moderate con- 
dition, i.e., driven by a weak force, in which the drift 
motion is submerged in the intense thermal diffusion and 
thus difficult to detect.[!4.19-24 This difficulty in deter- 
mining the drift velocity severely retards further explo- 
ration of the transport properties of the system working 
in an NESS. Besides, compared with the well-established 
equilibrium statistical mechanics theory, our understand- 
ing of non-equilibrium states is still very primitive. There- 
fore, until now a better understanding on the properties 


of systems in an NESS is very challenging from both the- 
oretical and practical points of view. 

Transport process is essential in the study of an NESS. 
Usually, particles in solids and liquids vibrate around their 
local equilibrium positions, and occasionally hop from 
one equilibrium position to another.!??] Therefore, trans- 
port properties such as mobility, conductivity, and self- 
diffusion can be simplified as a sequence of particle hops 
between effective potential wells.?°-?° The motion of a 
particle in a one-dimensional (1-D) periodic potential well 
is a simple but important model in studying various trans- 
port behaviors in an NESS. In particular, a large number 
of phenomena in physics and chemistry, such as surface 
diffusion, superionic conductors, rotation of molecules in 
solids, and charge carrier transport, can be understood 
based on this simple model.!?7] Compared with the in- 
tensive studies on the transport under a strong external 
force,!28—35] the properties and the underlying microscopic 
dynamics of the transport behavior in the weak force 
regime are still not fully understood, heavily attributed 
to the difficulty of determining the drift velocity in this 
regime. 

In the present paper, we employ a numeric method 
to accurately calculate the drift velocity of a particle 
driven by a weak external force in a 1-D periodic potential 
by non-equilibrium molecular dynamics (MD) simulation. 
Mobility and self-diffusion coefficient of the system can 
then be determined according to the obtained drift veloc- 
ity. We found that mobility is independent of external 
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force but increases with temperature. Moreover, it devi- 
ates from the Arrhenius law and follows a modified Ar- 
rhenius equation with its activation energy much smaller 
than the real energy barrier. Similarly, self-diffusion co- 
efficient also deviates from the Arrhenius behavior and 
along with mobility obeys the Einstein relation. A micro- 
scopic hopping mechanism was proposed to explain the 
non-Arrhenius behavior. Moreover, the position distri- 
bution of the particle in the potential well was found to 
deviate from the Boltzmann equation in an NESS. The 
non-Boltzmann behavior may be attributed to the ther- 
mostat introducing an effective “viscous” drag opposite to 
the drift direction of the particle, which also explains the 
fact that the activation energies for mobility and diffu- 
sion are much smaller than the real energy barrier of the 
potential well. 


2 Computational Procedures 


2.1 Drift Velocity, Mobility, and Self-Diffusion 


Drift velocity is an essential quantity for characteriz- 
ing transport properties. However, since the weak drift 
motion is typically submerged in the intense thermal 
diffusion, drift velocity that driven by a weak external 
force can hardly be detected in both experiment!®® and 
simulation.!!!:2°] For example, in a commonly used ionic 
liquid 1-ethyl-3-methylimidazolium tetrafluoroborate, the 
ion drift velocity is around 1077 m/s under a weak elec- 
tric field of 100 V/m, whereas the thermal velocity of ion 
can be as high as 190 m/s at 298 K.?4] This is the reason 
why most non-equilibrium simulations were carried out 
at a strong external force. There are three major meth- 
ods for calculating the drift velocity in an NESS: average 
velocity, mean-square displacement (MSD), and transient 
time-correlation function (TTCF). When the drift veloc- 
ity is comparable to or larger than thermal velocity under 
a strong force, it can be directly obtained from the aver- 
age velocity of all particles.!!1-37] Another way to calculate 
the drift velocity is comparing MSD with and without the 
external force. By using this approach, Murad!!?) calcu- 
lated the drift velocity of an NaCl solution under a strong 
electric field. This approach is only valid when the MSD 
is independent of an external force, but in many systems 
the force-dependent MSDs have been observed.!?!:38—39] 
Delhommelle, et al.!!9! employed the TTCF method to 
calculate the current density and conductivity of a molten 
salt under an external electric field. They found that the 
TTCF method can be used to analyze the response of a 
system to an arbitrary external force. However, in the 
TTCF approach an enormous number of non-equilibrium 
trajectories are required to perform data analysis, which 
limits the application of the TTCF method. 

In our previous work, a method has been successfully 
applied to calculate the drift velocities of ionic liquids and 
aqueous solutions under external electric fields.2!-38] Com- 
pared with the above approaches, this method only needs 
several trajectories and is able to obtain the drift velocity 


Vol. 62 


accurately at a moderate external force. Here we briefly 
describe this method as follows. The displacement of a 
particle under an external driving force in a time interval 
At at a time t; can be decomposed into a unidirectional 
drift term and a random diffusive term: 

Ax(ti) = ATarift (ti) + Avaitusion (ti) , (1) 
where Avarift (ti) and Awaiftusion (ti) are the displacements 
caused by drift motion and thermal diffusion, respectively. 
The two types of motions have totally different behav- 
iors: the drift motion shows a steady and unidirectional 
movement, whereas the thermal diffusion performs ran- 
dom walks and has a spectrum of white noise. For an 
infinite number of spontaneous displacements, we have 


N 

; 1 

Nim, N D Aaiffusion (ti) =0 ; (2) 
— 

where N is the number of time intervals. Therefore, the 

drift velocity can be approximately obtained from a finite 

length trajectory by the following expression 


PODES. 5 Atann (ti) _ 1S Arti) (3) 
PEEN. At N At ' 


i=1 i=1 
In an NESS the drift velocity remains constant, so we cal- 
culate a series of Ax with respect to various time intervals 
At and obtain viraj by linear fitting Ax against At. If we 
replace the long time average by the ensemble average over 
many independent short trajectories, we estimate the drift 
velocity by the equation, 

Varift © (Utraj) ; (4) 
-) denotes the ensemble average over all trajec- 


where (-- 
tories. 
The mobility u, characterizing the ability of ion move- 
ment in response to an external driving force f, can be 
calculated from the drift velocity by 
u = varitt/f - (5) 
In equilibrium, the self-diffusion coefficient D in one di- 
mension can be determined by using the Einstein expres- 
sion 
(Ag? (t)) = 2Dt, (6) 
where (Ax?(t)) is the MSD. This Einstein expression sug- 
gests that the MSD increases linearly with time for a dif- 
fusive motion. However, in an NESS the drift motion de- 
viates the MSD from a normal diffusive behavior. There- 
fore, the relative MSD that deducting the drift contri- 
bution from the overall MSD is defined to calculate the 
self-diffusion coefficient in an NESS as 
((Ax(t) — (Aa(¢)))?) = (Aa? (t)) — (Aa(t))? = 2Dt. (7) 
According to this equation, we can determine the self- 
diffusion coefficient in an NESS by linearly fitting the rel- 
ative MSD against time. 


2.2 Model 


The method described above can be easily applied to 
determining the drift velocity for various systems in MD 
simulation. In this work, a 1-D particle moving in a peri- 
odic potential was simulated to understand the behavior 
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of a system driven out of equilibrium by a weak exter- 
nal force. In our MD simulations, the particle moves in a 
bistable-potential model proposed by Sun"°l as 
V(x) = zf — 162”, (8) 
where x is the particle position and the periodicity is 
satisfied by applying the condition V(x + L) = V(x) 
with the periodic length L set to be 14. This model 
was used to study the biased sampling of non-equilibrium 
trajectories.4] In absence of the external force, the po- 
tential well is bistable with its minima separated by a 
low energy barrier Vo = 64 within the well and by a 
high energy barrier V; = 1681 between neighboring wells. 
Focusing on the weak force regime, the driving force f 
studied in this paper was chosen to be in the range from 
0 to 50 along the x+ direction and hence the potential 
well is effectively weakly tilted, as shown in Fig. 1. The 
temperature-dependent behavior was investigated at dif- 
ferent temperatures (T = 500, 750, 1000, 1250, and 1500), 
respectively, with the kinetic energy of the particle defined 
as r 
Eg = gel, (9) 


where the Boltzmann constant kz is scaled to be 1. 


Fig. 1 Tilted bistable effective potentials with different 
external forces in one period. 


2.3 Molecular Dynamics 


In the present paper, the velocity Verlet algorithm.!42] 
with a time step At = 0.001 was employed to integrate 
the Newton equation: 

dz dV (x) 

Wao a +f: (10) 

In order to detect the weak drift motion submerged in 
the thermal diffusion, an ensemble of 20 trajectories was 
sampled from 20 independent constant NVT MD simu- 
lations. Each simulation was carried out for 2 x 10°? MD 
steps with the particle position sampled every 1000 steps. 
The quantities studied in the present paper were averaged 
over all the 20 independent trajectories. The temperature 
was kept constant by using the Andersen thermostat.!43] 
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In the strong coupling limit, the particle frequently un- 
dergoes stochastic collisions and its new velocity is drawn 
from a Maxwell-Boltzmann distribution in every step. 


3 Results and Discussions 


3.1 Mobility 


Drift velocity is an important quantity for character- 
izing the transport process. However, as described above, 
it is very difficult to obtain the drift velocity under a weak 
driving force because the particle motion is usually dom- 
inated by the intense thermal diffusion. On the micro- 
scopic scale, if one tracks the movement of a particle un- 
der a weak external force in a short period of time, only 
isotropic random walks can be observed and the drift mo- 
tion is not detectable. Only when the observation time is 
sufficiently long can the hidden drift motion accumulate 
to be macroscopically measurable. Some typical trajecto- 
ries at T = 1000 under different external driving forces 
are shown in the inset of Fig. 2. It can be seen that the 
system exhibits a random motion caused by thermal dif- 
fusion and no noticeable unidirectional drift. However, 
no matter how weak the drift motion is, the microscopic 
structure and dynamics should have some signatures that 
are responsible to the macroscopic transport behavior. 


— f=0 4 
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f=40 ] 
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uni 4 
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Fig. 2 Average displacements as a function of time at 
T = 1000 under different external forces. Inset: Typical 
trajectories at T = 1000 under different external forces. 


In order to reveal the drift behavior submerged in 
the intense thermal diffusion, we applied the method de- 
scribed in Subsec. 2.1 to calculate the drift velocity in our 
simulations. Figure 2 shows the average displacements at 
different external forces. Compared with the real trajec- 
tories that show random walks, the average displacements 
increase linearly with time, showing a typical drift behav- 


ior. 
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Fig. 3 (a) Drift velocities and (b) mobilities at different 
temperatures and external forces. The solid lines in (a) 
depict the linear fits to the data. 


Figure 3(a) shows the drift velocities versus the exter- 
nal force at different temperatures and the solid lines are 
the linear fits to the drift velocities. By using the method 
described above, we managed to calculate the drift veloci- 
ties even at very small external forces (f L/Vı = 0.017). It 
can be seen that the drift velocity increases linearly with 
the external force, suggesting that the system remains in 
a linear-response regime. The drift velocity also increases 
with temperature at a given external force, but this trend 
becomes weaker as the temperature increases. The mobil- 
ity was then calculated from the drift velocity according 
to Eq. (5). Figure 3(b) shows the mobility against 1/T 
in a semi-logarithmic plot. As expected, mobility is inde- 
pendent of the external force in the linear-response regime 
and apparently increases with temperature. However, the 
mobility deviates from the Arrhenius law, 

u = u exp(— E; /ksT), (11) 
where u is the mobility, ug is the prefactor, E’, is 
the activation energy for mobility, ks is the Boltz- 
mann constant, and T is the temperature. The de- 
viation from the Arrhenius law has been observed in 
many systems. 
bility, conductivity, and diffusion coefficient, that devi- 
ate from the Arrhenius law can usually be described 


The transport properties, such as mo- 
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by the Vogel-Fulcher-Tamman (VFT) equation,[44-46] 
X = XoT7exp(—B/(T — To)) or a modified Arrhenius 
equation,47—5°] X = XT? exp(—E/ksT), where X rep- 
resents transport properties. In our system, we found that 
the mobility can be best fitted by the following equation, 


Eu 
u = wT? exp ( — a . 


The fitting results are listed in Table 1. The mean acti- 


(12) 


vation energy for mobility is 1052, much lower than the 
real energy barrier height of 1681. In the following sec- 
tions, we will show that the low activation energy and the 
non-Arrhenius behavior of mobility can be understood by 
a microscopic hopping mechanism. 


Table 1 The fitting parameters of mobility by Eq. (12). 


f Eu uo Ra 
0 = = à 

10 1077 0.12 0.9995 
20 1052 0.12 0.9996 
30 1049 0.12 0.9998 
40 1039 0.12 0.9990 
50 1041 0.12 0.9990 


a RÊ? is the coefficient of determination of the fit. 


3.2 Self-Diffusion 


In the linear-response regime, if the particles in a 
many-body system are uncorrelated, the mobility and self- 


diffusion coefficient satisfy the Einstein relation 
D=ukszT. (13) 


The self-diffusion coefficient was calculated according to 
Eq. (7). 
and self-diffusion coefficient at different temperatures and 


Figure 4 shows the relation between mobility 
external forces. It can be seen that the mobilities and 
self-diffusion coefficients indeed obey the Einstein relation 
very well within the statistical errors. 

According to Eqs. (12) and (13), the self-diffusion co- 
efficient also deviates from the Arrhenius law and can be 
well fitted by a modified version of the Arrhenius equation, 

Ep 
ET) 
where Do is the prefactor and Ep is the activation energy 


D = DoT"? exp ( = (14) 


for self-diffusion. Figure 5 shows the self-diffusion coeffi- 
cients against 1/T on a semi-logarithmic scale. It can be 
seen that the diffusion coefficient indeed can be described 
by Eq. (14) very well. The fitting parameters are given 
in Table 2. 
be 1048, very close to the activation energy for mobility 


The mean activation energy was found to 


and also much lower than the real energy barrier height 
of 1681. 
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Table 2 Fitting parameters of self-diffusion coefficient 
by Eq. (14). 


F Ep Do R? 

0 1021 0.12 0.9999 
10 1064 0.13 0.9992 
20 1068 0.13 0.9999 
30 1036 0.12 0.9991 
40 1013 0.12 0.9998 
50 1036 0.12 0.9995 

0.8 T T 7 T T T T 


— Einstein relation 

+ T=500 

+ T=750 

+ T=1000 
T=1250 
T=1500 


ukgT 


“0.0 0.2 0.4 0.6 0.8 


D 
Fig. 4 Einstein relation connecting the mobility and 
self-diffusion coefficient. For each temperature, different 
data points represent the data at different external forces. 
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Fig. 5 Self-diffusion coefficients at different tempera- 
tures and external forces. 


Up to now, we have successfully determined the mobil- 
ity and self-diffusion coefficient of the system driven out 
of equilibrium by a weak external force. We found that 
both mobility and self-diffusion show non-Arrhenius be- 
havior, but still obey the Einstein relation. To further 
understand the non-Arrhenius behavior, the distribution 
of particle positions (DPP) and hopping process are stud- 
ied in detail in the following sections. 


3.3 Distribution of Particle Positions in the 
Potential Well 


In equilibrium states, at a given temperature the DPP 
in a potential well obeys the Boltzmann distribution 


1 E 
p= beo (i 
zP URT” 


(15) 
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where kps is the Boltzmann constant, T is the tempera- 
ture, E is the potential energy, and Z is the normalization 
constant. When an external force is applied, the particle 
moves along the force direction with a drift velocity and 
the system reaches an NESS. Unfortunately, no general 
theories for the DPP in an NESS are currently available. 
In the present model, we calculated the DPPs at different 
external forces and temperatures. Figures 6(a) and 6(b) 
show the equilibrium DPPs at different temperatures and 
the non-equilibrium DPPs under various external forces 
at T = 500, respectively. The DPPs at other tempera- 
tures and external forces are not shown since they have 
no qualitative differences from those plotted in Fig. 6. 
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Fig. 6 Distributions of particle positions in the poten- 
tial well at (a) different temperatures at f = 0 and (b) 
different external forces at T = 500. 


As shown in Fig. 6(b), the external force obviously 
tilts the distributions in an NESS. We fitted the DPPs 
phenomenologically and found that the DPPs can be well 
described by a Boltzmann-type distribution function, 

1 p=(V(e) = f*2) 
P(x) = Z exp | kT” 

where P(x) is the DPP probability, V(x) is the poten- 
tial energy, Z is the normalization constant, f* and T* 
are phenomenological parameters. The fitting results at 
T = 500 are listed in Table 3. (The DPPs at other tem- 
peratures and forces can also be well fitted by Eq. (16) 
and the fitting results are not shown.) 


(16) 
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Table 3 The fitting parameters of particle position 
distributions by Eq. (16) at T = 500. 


f f* fr T Z R? 

0 —0.48 0.48 501.6 11.34 0.9996 
10 5.97 4.03 500.9 11.34 0.9996 
20 12.02 7.98 502.9 11.37 0.9995 
30 18.61 11.39 503.3 11.42 0.9995 
40 24.50 15.50 504.9 11.49 0.9994 
50 31.06 18.94 506.2 11.57 0.9992 


It can be seen that the phenomenological tempera- 
ture T* is approximately equal to the real temperature 
T. However, the phenomenological driving force f* is 
much smaller than the applied external force, indicating 
a deviation from the equilibrium behavior. The difference 
fr = f — f* is also listed in Table 3. We found that fr 
increases roughly linearly with the external force and sur- 
prisingly, strongly correlated with the drift velocity. Fig- 
ure 7 shows the relation between fr and drift velocity. It 
can be clearly seen that fr also increases linearly with the 
drift velocity and the slope changes slightly with temper- 
ature. All these features suggest that fy actually acts as 
an effective “viscous” force opposite to the drift direction. 


T x T T T j T 7 T 7 T J 
+ T=500 

40 F + T=750 ® 7 
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Jails + T=1250 a | 

T=1500 : 

r g r 
fr 204 m 4 
s: 

10 F i £ 7 
OF É u 
0.000 0.005 0.010 0.015 0.020 0.025 

Udrift 


Fig. 7 Relation between the effective “viscous” force 
and drift velocity. For each temperature, different data 
points represent the data at different external forces. 


To quantify the “viscous” force, we introduced an ef- 
fective “viscosity” € and fitted fr as a function of drift 
velocity by the following equation, 

fr = varitt « (17) 
The fitting results are given in Table 4. Clearly, the ef- 
fective “viscosity” increases slightly with temperature, in- 
ferring that, in our model the Andersen thermostat ex- 
erts an effective “viscous” force on the system to balance 
the external driving force and thus stabilize the system. 
From a microscopic point of view, the Andersen thermo- 
stat dissipates the system energy through stochastic col- 
lisions of particles with the reservoir. In an equilibrium 
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state, the collision is isotropic — the probability of a par- 
ticle undergoing a collision from right and left sides are 
the same. However, when we apply the Andersen ther- 
mostat to an NESS, the isotropic collision tends to elim- 
inate the anisotropic drift motion, since the symmetric 
Maxwell—Boltzmann distribution is assumed in the An- 
dersen thermostat, which brings in an effective “viscous” 
drag opposite to the drift direction. Here we should dis- 
tinguish the effective “viscous” force fr from the physical 
viscous force. The former effectively responds to the drift 
velocity in an NESS, whereas the latter is in response to 
the spontaneous velocity in a real viscous system. 


Table 4 Fitting parameters of effective “viscous” force 
by Eq. (17). 


T 500 750 1000 1250 1500 
E 1618.6 1708.6 1770.6 1817.4 1846.6 
R? 0.9997 0.9999 0.9999 0.9999 0.9999 


Based on the results above, we propose that the ther- 
mostat in an NESS seems to play an essential role in de- 
termining the non-equilibrium DPP. By introducing an 
effective “viscous” force in response to the unidirectional 
drift motion, the DPP in an NESS can be well described 
by the modified Boltzmann distribution function. 


3.4 Hopping Dynamics 

An NESS is always related to the transport process 
with a constant mass, charge, or energy flux. At the 
atomistic scale, the transport mechanism can be described 
by a sequence of incoherent hops of particles, which have 
been observed in many systems, such as liquid crystal,[51] 
ionic liquid,©?! glass forming liquid,!53] membrane,54 and 
polymer.®°! Therefore, it is of great importance to un- 
derstand the hopping mechanism behind the drift behav- 
ior. Here we analyze the hopping dynamics in our sys- 
tem by using the intermittent time correlation functions 
(TCFs):[56-57] 


(18) 


where the population variable p(t) is unity when the par- 
ticle stays in a particular potential well at time t and zero 
otherwise. The bracket (---) indicates an average over all 
starting time. By definition, C(t) describes the structural 
relaxation of a particle stays in the same potential well at 
time t as at time 0 after a sequence of hops. Because a 
particle moving away with a faster drift speed has a lower 
probability to travel back to its original position, the hop- 
ping dynamics C(t) can be used to study the mobility of 
the particle. Figure 8(a) shows C(t) in equilibrium at dif- 
ferent temperatures. It can be seen that C(t) decays much 
faster at a higher temperature, but this trend becomes 
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weaker as the temperature increases, consistent with the 
behavior of drift velocity (see Fig. 3(a)). 
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4000 
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Fig. 8 Intermittent time correlation function of hopping 
dynamics C(t) at (a) f = 0 and (b) T = 500. 


Figure 8(b) displays the non-equilibrium C(t) at T = 
500 under various external forces (C(t) at other temper- 
atures and external forces show no qualitative differences 
and are therefore not shown). Because a faster drift mo- 
tion can effectively decrease the probability of the particle 
moving back to its original well, C(t) decays faster under 
a stronger external driving force. Moreover, comparing 
the equilibrium C(t) with non-equilibrium C(t), one can 
find that the temperature and external force influence the 
hopping dynamics at two different time scales. As shown 
in Figs. 8(a) and 8(b), C(t) begins to change with tem- 
perature at a very short time scale, whereas the external 
force manages to affect the hopping dynamics at a longer 
time scale. This can be understood by the fact that tem- 
perature changes particle hopping rate immediately and 
hence influences C(t) at a short time scale. In contrast, 
since the external force is relatively weak, it can only break 
the symmetry of hops: the particle has a slightly higher 
probability to hop out of the well toward the force direc- 
tion and a slightly lower probability to hop backward. At 
a short time scale, the asymmetry of hops is negligible 
and as a result, the particle exhibits a random diffusive 
behavior. At a long time scale, the asymmetric hops ac- 
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cumulate to be detectable, leading to the decrease of C(t) 
and resulting in a unidirectional drift motion. 


3.5 Hopping Rate 


In the previous section, the force and temperature de- 
pendences of hopping dynamics were discussed by using 
the intermittent TCF. However, the TCF is insensitive 
to the hopping direction which is essential to the unidi- 
rectional drift motion. Therefore, the directional hopping 
rate defined as the number of hops toward each direction 
per unit time was calculated to further investigate the mi- 
croscopic dynamics of the system. In an equilibrium state, 
the particle performs random but symmetric hops in the 
potential well. When an external force is applied, the 
symmetry is broken and the hopping rates toward differ- 
ent directions become asymmetric, resulting in a net drift 
current. In our model, after being trapped for a period 
of time, the particle may escape from a potential well for- 
ward with a hopping rate kų or backward with a rate k_, 
and k = ki + k_ is the total hopping rate escaping out of 
the well. By definition, the drift velocity is connected with 
the hopping rate difference by the following equation, 


Varitt = (ky — k_)L*, (19) 


where L* is the mean hopping length. 


0.025 T T T T T T T 
| — vant = (k+ — k-)L* í 
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0.000 


0.0015 


0.0005 0.0010 


ka — k- 
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Fig. 9 Relation between the drift velocity and the hop- 
ping rate difference. For each temperature, different data 
points represent the data at different external forces. 


Figure 9 plots the relation between the drift velocity 
and the hopping rate difference. It can be seen that, at 
different forces and temperatures, the drift velocities and 
the hopping rate differences follow Eq. (19) very well. The 


mean hopping length L* was determined to be 14.1+0.8, 
which equals the well length within the statistical error. 
Therefore, according to Eq. (19), the transport property of 
the system can be further understood by the microscopic 
hopping process. 
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1000/T 
Fig. 10 (a) Hopping rate without the external force and 
(b) reduced hopping rate differences at various external 
forces as a function of temperature. The solid lines in (a) 
and (b) depict the best fits to the data by Eqs. (20) and 
(24), respectively. 


The hopping rate without the external force is shown 
in Fig. 10(a). 
with kip = k-o = ko, and the hopping rate increases 


In equilibrium the hops are symmetric 


significantly with temperature, consistent with the hop- 
ping dynamics discussed above. For a simple system, 
it is usually accepted that the hopping rate obeys the 
Arrhenius law,?2?:25.58-59] whereas a number of systems, 
such as glass,/48:60-61] ionic conductor,62-®! polymer 
electrolyte,!44-45:49] ionic liquid,!4°° and glass-forming 


d,[67—68] were found to show non-Arrhenius behavior. 


liqui 
In our model, we found that the equilibrium hopping rate 
also deviates from the Arrhenius law and can be well de- 
scribed (with R? = 0.99999) by a modified version of the 
Arrhenius equation similar to the one for self-diffusion, 
ko = AT exp ( = +) (20) 
where A is the prefactor and EF is the activation energy for 
hopping. The activation energy Ep was found to be 1528, 
very close to the real barrier V; = 1681. Here it is note- 
worthy that the real barrier V; is the maximum height be- 
tween the peak and valley of the potential well. However, 
since the particle has a probability to stay in any position 
of the well in a finite temperature (see Fig. 6), the activa- 
tion energy for a particle to hop out of the well should be 
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somehow lower than the maximum barrier height V;. Al- 
though the equilibrium hopping rate deviates from the Ar- 
rhenius law, the activation energy obtained from Eq. (20) 
turns out to be a good estimation of the energy barrier. 
Figure 10(b) shows the reduced hopping rate differ- 
ences against 1/T on a semi-logarithmic scale. The re- 
duced hopping rate difference is defined as the difference 
between the forward and backward hopping rates divided 
by its equilibrium value ko. If we assume that the external 
force changes the hopping rate by modifying the activa- 
tion energy, we can calculate the non-equilibrium hopping 
rate by modifying its equilibrium formula (Eq. (20)) by, 


k+ = AT? exp ( = —_— 


= too (95%) wh gh) 


where a is a parameter that describes the effect of exter- 
nal force on the activation energy. Therefore, the reduced 
hopping rate difference can be expressed as 
ky -k-  afh 
ko | kT? 
By fitting the simulation data, we found that a exponen- 
tially depends on temperature, 


Ea 
Q = Qo exp (7 ) é 
According to Eqs. (22) and (23), the reduced hopping rate 
difference can be well fitted by the following equation 
ky —-k-  aofL ( Ea ) 
= exp ; 
ko ks T ks T 
The fitting results are displayed in Fig. 10(b) and given in 
Table 5. Combining Eqs. (19), (20), and (24), we have 


(21) 


(22) 


(23) 


(24) 


k} —k_ 
varift = (k4 — k-)L* = kg 1" 
0 
2 

z a T-?/3 exp (- Fx — Ea) 

B B (25) 
Therefore, the mobility can be expressed as 
Varit  AaoL? -2/3 ( Ex — 2 
= = ——T —_——.— } . 2 

F k, exp ET (26) 


This equation explains why the mobility is independent 
of the external force, and gives the exact (Non-Arrhenius) 
temperature dependence as we have already obtained from 
Eq. (12). Comparing Eq. (12) with Eq. (26), we found that 
the prefactor (Aao L? /ks) equals uo and the activation en- 
ergy (Ek — Ea) is 1041, very close to Eu of 1052. This 
suggests that, by introducing a temperature dependent 
parameter a, the dependence of mobility on temperature 
and force can be well described by a microscopic hopping 
mechanism. Moreover, the exponential dependence of a 
on temperature explains the fact that the activation en- 
ergy E„ becomes much lower than the real energy barrier, 
considering that Ex is a good estimation of the real energy 
barrier. 


No. 4 
Table 5 The fitting parameters of reduced hopping 
rate differences by Eq. (24). 
F ao Ea R? 
0 z x z 
10 0.15 468 0.9999 
20 0.15 478 0.9999 
30 0.15 493 0.9999 
40 0.15 498 0.9999 
50 0.15 499 0.9999 


500 750 


1000 1250 1500 


T 


Fig. 11 Parameter a calculated by Eqs. (23) and (27) 
at different temperatures and external forces. 


The parameter a was usually chosen to be 1 in theoret- 
ical calculations, !??-6°! but in our model we found that it is 
actually smaller than 1 and becomes temperature depen- 
dent in an NESS. Figure 11 plots «a calculated by Eq. (23) 
in red lines. It can be seen that œ decreases with temper- 
ature and takes a value between 0.4 and 0.2. Moreover, a 
seems to be force independent. According to Eq. (21), af 
can be considered as an effective force that drives the par- 
ticle hopping. On the other hand, since we have described 
in Subsect 3.3 that the Andersen thermostat exerts an ef- 
fective “viscous” force fr on the particle, the net effective 
force experienced by the particle can also be expressed as 


af =f- fr. (27) 
According to Eq. (27), the parameter a was also calcu- 
lated from the effective “viscous” force fyr and shown in 
black lines in Fig. 11. It can be seen that the effective 
driving forces af obtained from DPP and hopping rate 
display similar force and temperature dependences. The 
quantitative difference may be attributed to the assump- 
tion in dealing with the effect of force on the hopping 
rate. For example, in Eq. (21) the external force changes 
the hopping rate by simply modifying the activation en- 
ergy with +a fL/2, but in fact the effect of the external 
force is expected to be more complex, since it depends 
on the DPP in the potential well. Moreover, the decrease 
of a with temperature can be naturally explained on the 
basis of the “viscous” force: the “viscosity” € and drift 
velocity vari, both increase with temperature, and as a 
result the effective driving force af decreases with tem- 
perature due to the increasing “viscous” force fr = Evaritt 
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(Eq. (17)). This again supports our previous conclusion 
that, in an NESS, the Andersen thermostat exerts an ef- 
fective “viscous” force on the system. As a result, the 
transport properties of the system can be understood as a 
sequence of non-Arrhenius hops under an external driving 
force and an effective “viscous” drag. 


4 Conclusion 

In conclusion, by non-equilibrium MD simulation, the 
drift velocity of a 1-D particle moving in a periodic poten- 
tial well under a weak external force has been successfully 
determined. The mobility and self-diffusion coefficient of 
the particle were also obtained by using the drift veloc- 
ity. We found that the mobility is independent of the 
external force and increases with temperature, and the 
mobility deviates from the Arrhenius law and follows a 
modified Arrhenius equation with the activation energy 
much smaller than the real energy barrier. Similarly, the 
self-diffusion coefficient also deviates from the Arrhenius 
behavior. Nevertheless, the mobility and the self-diffusion 
coefficient still obey the Einstein relation. 

To further understand the microscopic structure and 
dynamics in an NESS, the DPP and hopping dynamics 
in the potential well were calculated. We found that the 
non-equilibrium DPP can be well described by a modified 
Boltzmann distribution function if we introduce an effec- 
tive “viscous” force which probably originates from the 
energy dissipation by the Andersen thermostat. More- 
over, the microscopic hopping dynamics was analyzed by 
the intermittent TCF. In contrast to the temperature that 
changes the hopping rate at a short time scale, a weak ex- 
ternal force can only slightly break the symmetry of hop 
that accumulates to result in a unidirectional drift motion 
at a long time scale. We also propose a microscopic ex- 
planation to the non-Arrhenius behavior of mobility by 
a hopping mechanism. By introducing a temperature- 
dependent effective driving force, we successfully explain 
the fact that the activation energy E„ for mobility be- 
comes much lower than the real energy barrier. This again 
supports our result of the DPP, inferring that the Ander- 
sen thermostat exerts an effective “viscous” drag on the 
particle opposite to the drift direction. 

This work advances our understanding of the transport 
properties and processes in an NESS. Even though all the 
results in this work are based on a 1-D abstract model, the 
method of determining the drift velocity from diffusion- 
dominant trajectories and the microscopic hopping mech- 
anism may be easily extended to investigate other abstract 
models and real materials in an NESS. Since our work 
only studied a 1-D abstract model coupled to the Ander- 
sen thermostat, whether or not these non—Boltzmann and 
non-Arrhenius behaviors are general in other models and 
thermostats is an interesting topic for future work. 
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